Multihistogram Reweighting for Nonequilibrium Markov Processes Using Sequential 

Importance Sampling Methods 
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We present a multihistogram reweighting technique for nonequilibrium Markov Chains with dis- 
crete energies. The method generalizes the single histogram method of Yin et al. [Phys. Rev. E 
72, 036122 (2005)], making it possible to calculate the time evolution of observables at a posteriori 
chosen couplings based on a set of simulations performed at other couplings. In the same way as 
multihistogram reweighting in an equilibrium setting improves the practical reweighting range as 
well as use of available data compared to single histogram reweighting, the method generalizes the 
multihistogram advantages to nonequilibrium simulations. We demonstrate the procedure for the 
Ising model with Metropolis dynamics, but stress that the method is generally applicable to a range 
of models and Monte Carlo update schemes. 



I. INTRODUCTION 

In the last two decades the use of short time critical 
dynamics (STCD) as a way of investigating the critical 
properties of models in statistical mechanics has emerged 
as an interesting alternative to equilibrium Monte Carlo 
(EMC) simulations. The reason is twofold: First, STCD 
simulations may, given sufficiently large system sizes, ef- 
fectively avoid finite size effects, since the correlation 
length is still small in the short time regime. Second, and 
even more important, critical slowing down is avoided 
due to the system being far from equilibrium. See for 
instance Refs. 1 and 2 for some recent reviews on the 
STCD method. 

The EMC simulations have, on the other hand, long 
had the advantage of the powerful histogram reweight- 
ing methods of Ferrenberg and Swendsen, meaning 
that knowledge of the exact couplings of interest are 
not needed a priori the simulations: Single histogram 
reweighting [ ] makes it, in principle, possible to extract 
observables at arbitrary couplings from a simulation per- 
formed at a fixed coupling, although the practical range 
is limited to a small neighborhood of the simulation cou- 
pling. The practical reweighting range as well as the 
statistical quality of the observables can be improved by 
using multihistogram reweighting [4], where data from 
an arbitrary number of simulations at different couplings 
are combined in an optimal way. 

A first attempt to introduce reweighting techniques to 
STCD simulations was done by Lee and Okabe [ ] . They 
proposed what essentially is (Monte Carlo) time depen- 
dent single histogram reweighting, where the reweighting 
is done simultaneously with the simulation. This is a ma- 
jor drawback of the method, as all the desired couplings, 
simulation as well as reweighted, must be known a priori. 
In addition, obtaining an observable through reweighting 
is almost as computationally heavy as performing the 
simulation itself, and since the procedure scales linearly 
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with the number of reweighting couplings, the practical 
usefulness of the method is limited. 

Yin et al. proposed a novel method that improves 
the single histogram reweighting for STCD, but the 
method is restricted to models with a discrete energy 
spectrum [ ] : By keeping track of a time dependent his- 
togram of the energy changes at each step in the Markov 
chain, and sampling these each time an observable is sam- 
pled, it is possible to efficiently reweight observables to 
a posteriori chosen couplings. This is a significant im- 
provement, but the reweighting is still limited to the sin- 
gle histogram regime where only data from simulations at 
one coupling strength is used. In a real world situation 
one would typically perform several simulations at dif- 
ferent couplings when searching for the critical one, thus 
obtaining a lot of data that contains potentially useful in- 
formation. It would therefore be advantageous to develop 
a multihistogram reweighting technique for nonequilib- 
rium Markov processes. 

In this paper we propose such a technique for nonequi- 
librium reweighting of models with discrete energy spec- 
tra by extending the method of Yin et al. to time de- 
pendent multihistogram reweighting using the method 
of Ferrenberg and Swendsen. After the general deriva- 
tion we test the method on a Metropolis dynamics STCD 
simulation on the 2D Ising model and compare it with 
the "raw" (non-reweighted) as well as single histogram 
reweighted results. 

Note that although these methods are developed with 
STCD simulations in mind, they are generally applicable 
to all nonequilibrium Markov processes (and hence also 
equilibrium processes) . 



II. MULTIHISTOGRAM REWEIGHTING 

A. General Formulation 

Let Xt = (<7i, <T2, . . . , at) denote a Markov chain of field 
configurations a after t steps (often called "time" from 
now on), and let E t be the dynamical phase space, i.e. the 



set of all possible Markov chains of t steps, starting from 
an initial set of states So. Given a (usually small) set of 
coupling parameters f3, we may associate a probability 
weight to each Markov chain: 
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We now turn to the dynamic generalization of this 
method. Let X q> t = { x Qt t } be a set of N q Markov chains 
obtained from independent simulations at the same cou- 
pling f3 q , and let X t = { X q _ x \ q G {1, . . . , Q} } be the 
collection of all Markov chains from Q such sets, each 
(1) with their own associated coupling j3 q G {/3i, . . . , /3q }. 



At each time step we have a probability distribution over 
all possible reachable states, which we may capture in the 
dynamical partition function 



W) = £>(&**), 



(2) 



where x t under the summation sign is a shorthand nota- 
tion for x t G £4. In the same way we define the dynamical 
average of an observable O as 



(O^fflsZtW-^wfaxtMat) 



(3) 



Thus, at a fixed time t there is no difference between the 
dynamical formalism and the equilibrium formalism, ex- 
cept that the fundamental "states" in the dynamical case 
are entire Markov chains instead of single field configu- 
rations. With this in mind, methods from equilibrium 
simulations may be carried right over to dynamic simu- 
lations. 

For the sake of simplicity, we will from now on re- 
strict ourselves to the special case when the set of cou- 
pling parameters consists of only one coupling parameter 
(also named j3). Generalization to multiple parameters 
is straightforward. 

In the equilibrium formalism, Ferrenberg-Swendsen 
multihistogram reweighting goes as follows [4]: Let 
H(a) — Ho(a) + (3S(a) be a general Hamiltonian of a 
field a and the partition function Z(/3) — J2a exp(P(cr)). 
S is an operator (energy, magnetization, etc.) on the 
field. Given Q EMC simulations performed at couplings 
{ /3i, . . . , j3 q , . . . , Pq }, each with N q samples with auto- 
correlation time T q , the best approximation for the prob- 
ability weight of S at a (a posteriori chosen) coupling /3* 
is given by 



P(f3*,S) = 



E g (l + 2T,)- 1 Af g (S)exp(/rS) 
V N q (l + Zr,)- 1 exp(/3 9 S)Z(/y-i ' 



(4) 



where M q (S) is the histogram of the S- measurements of 
the g'th series, and 



Z{0)=J2 P ^'^- 



(5) 



P is found by solving Eqs. (4) and (5) iteratively. The 
error minimized average of an observable (operator on S) 
at f3* is then given by 



Replacing exp(/35 l ) — > w((3,x t ), S — ¥ x t (i.e. x t may be 
seen as the identity operator), P(/3, 5) — > W((3,x t ), and 
Z((3) -> Z t {fi) = E Xt ex t W((3,x t ) in Eq. (4), we get 
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j: q M q (x t )w([3*,x t ) 
E a N q w(/3 q ,x t )Z t (j3 q )- 



(7) 



r q = since the Markov chains are independent. M q (xt) 
is the histogram of Markov chains on the form xt for X qX , 
i.e. M q {x t ) — J2 y ex t &y,xf Now, if we treat all the ob- 
tained Markov chains as unique (even in the improbable 
case of two being identical) , there will be only one j/el ( 
fulfilling y = x t , namely Xt itself. Hence we must have 
that ^2 q M q (xt) = 1, and Eq. (7) simplifies to 



W(0*,x t ) 



w(j3*,x t ) 



Z q N q w(l3 q ,x t )Z t (f3 q )-i 
The error minimizing dynamical average is then 



(8) 



(0) t {P*) = Zti/3*)- 1 E 0(a t )W(f3*,x t ), (9) 

which is what we are interested in. To get any further 
with a practical use of Eqs. (8) and (9) we need to de- 
termine the Markov chain weights {w(f3,x t ) }■ 

B. The Markov Chain Weight 

The weight of a Markov chain [Eq. (1)] is per definition 
given by 



w(f3,x t ) = w(<7o) JJw()9,Ot|ot_i), 



(10) 



where w(cro) is the weight of the initial state and 
w(/3, at\(Tt-i) is the weight associated with a Markov step 
from o~t-i to at at coupling /?. Since each Markov chain 
in an STCD simulation may consist of a very large num- 
ber of steps (e.g. 10 10 for a 10 4 sweep simulation on a 
1000 x 1000 lattice), it is practically impossible to di- 
rectly store the weight information of all steps for later 
reweighting; in the general case, where the set of possible 
step weights is large (possibly infinite), one has no choice 
but to reweight "on the fly", as Lee and Okabe proposed. 
If, however, we restrict ourselves to models where the 
set of possible step weights is discrete, it is possible to 
find the the Markov chain weight at an a posteriori chosen 
coupling j3. The trick is to write Eq. (10) as [ ] 



(o)(/n ^(/n^E^p^s). 



(6) 



w(P,x t )=u J (o- )~[[u J ((3,ay" 



(11) 



w(/3, a) being the weight of a step of type a, and r a ,t G No 
the number of such steps in the Markov chain xt- w(/3, a) 
can be calculated for an arbitrary ft, given the knowledge 
of a and the dynamics used in the simulation. Thus, to 
know w((3,Xt) we only need to know the initial state <To 
and the histogram of r Qit -values. 

As an example, consider an STCD simulation of a spin 
model using canonical Metropolis dynamics at some cou- 
pling fi q . A trial update of a spin will result in an energy 
change AE a , which we assume to be member of a dis- 
crete set of all possible (one step) energy changes { AE }. 
There are three possible outcomes of a trial update in the 
Metropolis dynamics: (See Ref. 5 for mathematical de- 
tails.) 

• AE a < 0, in which case the trial is always accepted. 
Then u(f3, a) = 1. 

• AE a > and the update is accepted. Then 

uj((3,a) = cxp(-(3AE a ). 



• AE a > and the update is rejected. 
cj(ft,a) = l-exp(-(3AE a ). 



Then 



The explicit expression for the Markov chain weight, 
Eq. (11), becomes 

w(fi,xt)=u(a ) H [exp(-(3AE a )] r "- t 

a. 
AE a >0 
accept 

x [] [1 - exp(~f3AE a )] r ^ 

a 

AE a >0 

reject 

= w(a )exp(-/3A t £ a ) 

x H [l-exp(-ftAE a )] r ^\ (12) 

a 

AE a >0 

reject 



where 



A,E a 



E ' 

Q 

AE a >0 
accept 



tAE n 



E 

t 

A t E>0 
accept 



A t E 



(13) 



is the running sum of all positive acceptance energy 
changes. So, to be able to calculate the Markov chain 
Metropolis dynamics weight at an arbitrary /3, and thus 
to be able to multihistogram reweight an observable av- 
erage using Eq. (9), one has to keep track of A t E & as well 
as a histogram r Qit over the distribution of the rejected 
positive energy changes. For many discrete models this 
is a managebly small set of extra information to sample 
in addition to the standard obscrvables. 



III. TESTING THE METHOD ON THE 2D 
ISING MODEL 

We test the nonequilibrium multihistogram reweight- 
ing method on a Metropolis dynamics STCD simulation 
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FIG. 1. Plot of magnetization curves for the mul- 
tihistogram reweighted couplings (from below) /3* = 
0.439 25, 0.439 50, ... , 0.441 75. The reweighting is based on 
simulations at f3 q = 0.4400, 0.4405 and 0.4410. Line thick- 
ness corresponds to ± jackknife standard error. 



of a two dimensional ferrormagnetic Ising model on a 
square L x L = 64 x 64 lattice with periodic boundary 
conditions. The Hamiltonian is given by 



H\s\ = - 



E< 



H °J i 



(14) 



where i and j are lattice indices and Sj G {—1,1}. Note 
that {AE} = {-4,-2,0,2,4}, which means that we 
just need to sample three extra quantities in addition to 
the observable(s), in this case the dynamical average of 
the magnetization per site, 



rm^L^Y, 



St,i- 



(15) 



To make the test simple, we consider only simulations 
from a perfectly ordered state, So = Co = { $i = 1 | Vf }, 
and so we may choose w(ao) = 1- 

Figure 1 shows a fan of reweighted magnetization 
curves as a function of Monte Carlo time, obtained from 
three simulations performed at j3i — 0.4400, /3 2 = 0.4405, 
and /?3 = 0.4410. (This is close to criticality, (3 C = 
ln(l + \fl)l2 = 0.440 068 7. . .) N q = 20 000 Vq. Errors 
are obtained by the Jackknife method. Notice how the er- 
ror grows as the reweighted coupling deviates more from 
the simulated couplings and the overlap in energy his- 
tograms worsens, as is the case for all statistical reweight- 
ing techniques. 

We compare the multihistogram reweighting with raw 
data and single histogram reweighting. First we perform 
single histogram reweighting of the q = 1 simulation to 
/3* = j3 2 and from simulation 3 to (3* = (3 2 , then multihis- 
togram reweighting to f3* — p 2 using simulation 1 and 2, 
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and finally multihistogram reweighting to (3* = 02 using 
the entire set of available data, q = 1, 2, 3. The results are 
compared with the non-reweighted dynamical average of 
the raw data of simulation 2 in Fig. 2. Notice that while 
the single histogram reweighting always performs worse 
(in terms larger error bars) than the raw data, the mul- 
tihistogram reweighting may even outperform the raw 
data, given that, like in this case, more simulations with 
overlapping energy histograms exist. 

It should be noted that the computational cost of mul- 
tihistogram reweighting is of the same order of magnitude 
as single histogram reweighting - and negligible com- 
pared to the simulations. 

IV. CONCLUSION 



To summarize, we have generalized the single his- 
togram nonequilibrium Markov Chain reweighting tech- 
nique of Yin et al. [6] to a multihistogram framework 
based on the equilibrium method of Ferrenberg and 
Swendsen. In this manner we can take advantage of 
several noncquilibirum simulations performed at differ- 
ent (but close) couplings, producing dynamical averages 
with equal or smaller errors than those obtained by av- 
erages based on the non-reweighted or single histogram 
reweighted datasets. 



FIG. 2. Comparison between non-reweighted magnetizations 
curves (dashed gray lines) at /3 = 0.4405 and reweighted 
curves (thick black lines), reweighted to /3* = 0.4405. (a) 
shows the single histogram reweighting from a simulation at 
fiq = 0.4400. (b) shows the single histogram reweighting from 
a simulation at f3 q = 0.4410. (c) shows a multihistogram 
reweighting based on both these simulations, (d) shows a 
multihistogram reweighting using all data available, i.e. the 
simulations as fi q — 0.4400, 0.4405 and 0.4410. Error bars are 
jackknife standard errors. 
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